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Abstract 

We perform a joint analysis of the cosmic microwave background (CMB) and Galactic emission 
from the WMAP 7-year temperature data. Using the Commander code, based on Gibbs sampling, 
we simultaneously derive the CMB and Galactic components on scales larger than 1° with improved 
sensitivity over previous work. We conduct a detailed study of the low-frequency Galactic foreground, 
focussing on the "microwave haze" emission around the Galactic center. We demonstrate improved 
performance in quantifying the diffuse Galactic emission when including Haslam 408MHz data and 
when jointly modeling the spinning and thermal dust emission. We examine whether the hypothetical 
Galactic haze can be explained by a spatial variation of the synchrotron spectral index, and find that 
the excess of emission around the Galactic center is stable with respect to variations of the foreground 
model. Our results demonstrate that the new Galactic foreground component - the microwave haze - 
is indeed present. 

Subject headings: cosmology: observations, cosmic background radiation, diffuse radiation, methods: 
data analysis, numerical, statistical, Galaxy: center 



1. INTRODUCTION 

While data from t he Wilkinson Microwave Anisotropy 
Probe (WMAP) (see lJarosik et al.ll201tt iKomatsu etaLl 
120101 and references therein) has enabled unprecedented 
advances in the understanding of cosmology over the past 
decade, it has also opened a unique window into the fun- 
damental physical processes of the interstellar medium 
(ISM) . The choice of observing bands for WMAP insured 
that multiple emission mechanisms would be observed 
across the frequency coverage. In particular, there are at 
least three distinct physical processes at low frequencies 
(23-33-41 GHz),: free-free, synchrotron, and anomalous 
microwave emission (AME), which falls with frequency 
and is highly correlated with 100 (im thermal dust emis- 
sion. At high frequencies (61-94 GHz), Galactic emis- 
sion is completely dominated by thermal dust emission 
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(|Gold et alJ[20Tl . 

Free-free emission (or thermal brcmsstrahlung) is gen- 
erated by scattering of ionized electrons off the pro- 
ton nuclei in hot (~ 5000 K) gas, and has a bright- 
ness temperature which scales as T oc v through 
the WMAP bands, where v represents frequency. The 
bulk of the synchrotron emission observed by WMAP is 
seen t o closely follow a power-law T oc v~ 3 (jKogut et al.1 
120071 iDobleri 1201 lh . Lastly, spinning dust (whose 
presence has b een observed in small, dusty clouds 
ijCasassus et al.H2008t fPlanck Collaboration et al.ll2011l) 
as wel l as hotter diffuse regions ([Dobler fe Finkbeinerl 
I2008bl iDobler et aD 12009ft ). is the likely cause of the 
anomalous microwave emission. Small dust grains with 
non-zero dipole moments are spun up by a variety of 
mechanisms such as ion collisions, plasma density fluctu- 
ations, photon fields, etc., and produce spinning dipole 
radiation (see lEricksor] 1 19571 iDraine fc Lazarianl 119981 
for the original theoretical realization of this idea). 

Using simple templat e regre ssion t echniques (s e e Sec - 
tion EJ, iBennett et aD (|2003h and IFinkbeinerl ()2004l ) 
showed that the Galactic emissions are highly spatially 
correlated with maps at other frequencies. Free-free 
emission is morphologically correlated with Ha recom- 
bination line emission, synchrotron with low frequency 
radio emission (e.g., at 408 MHz), and spinning (and 
thermal) dust with total dust column density (e.g., 
ISchlegel et al.l (119981 ) evalua ted using models for the ther- 
mal emission to 94 GHz by [Finkbeiner et al.lll999 f). Af- 
ter removing emi ssion correlated with these templates, 
IFinkbeinerl (j2004l ) found that there was an excess sig- 
nal centered on the Galactic center (GC) and extend- 
ing out roughly ~ 30 degr ees. A more detailed study 
of this Galactic "haze" by (jDobler fc Finkbeinerl l2008a| 
(and more recently with the 7-vear data bv lDoblerll2011l i 
showed that its spectrum (T oc v& with (3 ~ —2.5) was 
too soft to be free-free emission and too hard to be syn- 
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chrotron emission associated with acceleration by super 
novae (SN) shocks (after taking into account cosmic-ray 
diffusion). The origin of this residual component remains 
a mystery and has been the object of intense theoreti 

^ ^^^^^^^^^^ ^ ^^^^^^^^^ 

cal scrutiny (e.g.. iHooper et al.l 120071: ICholis et al.ll200 
Biermann et al] 20101 iDobler et al.ll201lt iCrocker et al 
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However, the subject of the haze has not been with- 
out controversy, most notably due to t he claim by the 
WMA P team (as well as others, e.g., iDickinson et al.l 
I2009T ) that the haze is no t dete cted in their analyses. 
In particular, IGold et ail (|2011[ ) claimed lack of evi- 
dence for the haze based on WMAP polarization data, 
though it is un clear that a p olarization signal would be 
detectable (see IDobler] 1201 ID . A comprehensive study 
of compon ent separation (and CMB c l eaning) was per- 
forme d by lEriksen et all (|2006l . I2007bl) : IDickinson et al.l 
(|2009f ) utilizing Bayesian inference of foreground ampli- 
tudes and spectra via Gibbs sampling. These studies 
have never claimed a significant excess of emission to- 
wards the GC. Nevertheless, with the release of gamma- 
ray data from the Fermi Gamma-Ray Space Telescope 
in 2009, a corresponding feat ure at high energies was 
discovered (jDobler et al.ll2010l ). likely generated by the 
haze electron inverse Compton (IC) scattering starlight, 
infrared, and CMB photons up to Fermi energies. 

Wi th few exceptions (see for example iBottino et al.l 
2010;), foreground studies using template fitting pre- 
subtract a CMB estimate which imprints a bias in 
the foreground spectra (see IDobler fc Finkbeiner|[2008al) 
since no CMB estimate is completely clean of fore- 
grounds. In this work, we attempt to minimize this effect 
by solving jointly for the CMB map and angular power 
spectrum and Galactic emission parameters, within a 
Bayesian framework. We fit foreground parameters in 
every pixel, thus taking into account spatial variation of 
foreground spectra, and we test the stability of the so- 
lution against several Galactic emission models. Seeking 
greater independence from correlation with external tem- 
plates, we instead treat other data sets as input channel 
maps and process them through th e Gibbs sampler: to 
this end we add Haslam at 408 MHz (jHaslam et al.lll982t) 
to the WMAP data set. 

In Sj2]we describe the data and our foreground model, 
while in fj3] we describe our its application in our 
Commander analysis. Our approach enables us to improve 
upon previous studies. In 2] we discuss our results and 
refine the methodology by performing a joint analysis of 
the WMAP 7-year data set and Haslam at 408 MHz. 
Lastly in £[S]wc summarize and draw our conclusions. 

2. DATASETS AND FOREGROUND MODEL 

Our approac h uses the Gibbs samp ling algorithm 
introd uced by Uewell et all (120041) and IWandelt et al.l 
(120041) . and further d e veloped by Eriksen et all (120041); 
lO'Dwver etHI (l200l : l Eriksen et al.l (l2007bD; iChu et al.l 



(120051) : Uewell et al 



lLarson et all (|2007D 



(20091): IRudiord et al.l (120091) : 

The method has been 



merically implemented in a computer program called 
Commander, and successfully applied to previous re- 
leases o f the WMAP da t a as re ported in lEriksen et al 



2008al) . lEriksen et all (|2007aD . and IDickinson eta! 
2009). Commander is a maximum likelihood method 



which generates samples from the joint posterior density 



for the CMB map and angular power spectrum, as well as 
foreground components for the chosen sky model. A de- 
tailed description of the algor ithm and its validatio n on 
simulated data is provided by lEriksen et al.1 (|2007bl and 
references therein). For the interested reader, we sum- 
marize the technical aspects of the sampling procedure 
in Appendix [5] 

The main advantage of the approach is full charac- 
terization of the (posterior) probability distribution of 
the parameter space spanned by the adopted sky model. 
Moreover, we can evaluate the goodness-of-fit of each 
sample, i.e., of a given set of model parameters. Any 
parametric model can be encoded, leaving the freedom 
to vary parameters pixel-by-pixel as well as to fit tem- 
plates. 

In the following we describe the data used in our anal- 
ysis and their processing. 

2.1. WMAP 7-year dataset 

Our aim is to perform a compre hensive study of th e 
WMAP 7- year temperature data ([Jarosik et al.l 120101) , 
focusing not only on estimating the CMB signal and cor- 
responding power spectrum, but also on characterizing 
the properties of the foreground emission, with particu- 
lar attention to evidence for the haze signal. 

The WMAP data set comprises data from ten Differ- 
ential Assemblies (DAs) covering five frequencies from 23 
to 94 GHz. The WMAP team provided maps centered at 
frequencies of 23, 30, 40, 60 and 90 GHz [K,Ka,Q,V,W 
bands] resulting from the coadition of the DAs. We 
smooth the frequency maps to an effective 60 arcminutc 
resolution by deconvolving the maps with the appropri- 
ately coadded transfer functions provided by the WMAP 
team and convolving with a Gaussian beam of 1 degree 
FWHM. We then downgraded the sky maps to a work- 
ing resolution of N S id c = 128 in the HEALPix 12 scheme 
(|G6rski et al.ll2005l ). The choice of angular resolution, 
and subsequently of the pixelization, was dictated by 
both the angular resolution of the available foreground 
templates and by the need to resolve the first acoustic 
peak of the CMB power spectrum. This represents a 
novel element of our analysis compared to p revious work 
(|Eriksen et al.ll2008afc IDickinson et al.ll2009l) . 

In a further break with earlier analyses, we do 
not add uniform wh ite noise to regularize the maps 
(jEriksen et al.ll2007a| ). but instead add a noise compo- 
nent proportional to the actual noise variance in the 
smoothed, processed maps. To do so, we compute the 
RMS noise per pixel of the maps at the lower resolu- 
tion via a Monte Carlo approach. We generate 1,000 
non-uniform white noise realizations for each channel fol- 
lowing the prescription given by the WMAP team. In 
practice, we draw random Gaussian noise maps for each 
frequency band with zero mean and a variance given by 
crl/N ohs at the native WMAP resolution N sidc = 512, 
where N \, s represents the number of observations in a 
given pixel. We then smooth and downgrade the noise 
maps as above and finally recompute the variance per 
pixel averaged over 1,000 simulations. Noise is then 
added to each frequency based on the computed noise 
variance, and chosen so that the signal-to-noise ratio of 
the masked smoothed map at i ~ 2iV si d c is of order 

12 http://hcalpix.jpl.nasa.gov/ 
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unity. The scaling factors applied arc [16,12,12,6,6] for 
the [K,Ka,Q,V,W] channels respectively. One of the ad- 
vantages of this approach is that it preserves the noise 
structure imposed by the scanning strategy and describes 
at least the diagonal part of the instrumental noise, which 
has been correlated by the smoothing procedure. 

2.2. Foregrounds 

The diffuse Galactic emission consists of three contri- 
butions from well understood foregrounds - synchrotron 
radiation from cosmic ray electrons losing energy in the 
Galactic magnetic field, free-free emission in the diffuse 
ionised medium, and radiation from dust grains heated 
by the interstellar radiation field. In addition, there is 
a strong contribution at frequencies in the range 20 - 
70 GHz referred to as anomalous microwave emission 
(AME) that is strongly correlated with the thermal dust 
emission and which has been explained by rotational ex- 
citation of small grai ns - the so-called spinning dust emis- 
sion (see for example IHoang et al.ll2010l 120111 for recent 
studies on the subject and references therein) . In the fre- 
quency range spanned by WMAP, the synchrotron emis- 
sion is reasonably well described by a power law bright- 
ness temperature emissivity with spectral index j3 ~ —3. 
The free-free emission also follows an approximate power 
law, well described by a spectral index a ~ —2.15. The 
thermal dust component is usually modeled as a gray 
body with emissivity of the form T{y) cx v e B(y, Td) with 
typical values for the parameters given by e ~ 1.6 — 2.0 
and Td ~ 18K. Since the highest frequency WMAP band 
W has a nominal central frequency of 94 GHz, the over- 
all thermal dust contribution is small and can be ap- 
proximated by a simpler power law model with a spec- 
tral index e ~ 1.7. The AME is not completely well- 
characterised as yet, but falls rapidly below 20 GHz. Fits 
to the high latitude sky suggest that it may be reasonably 
well approximated by a power-law emissivity over the 
WMAP range of wavelengths, although this may reflect 
the combination of multiple spinning dust populations in 
different physical conditions along a given line-of-sight. 

The total foreground intensity observed by the v chan- 
nel in a given direction p can then be summarized as 
follows: 

D v<d (p)+[ — ) A syncb {p) 

d—x,y,z 

( V \ Q(P) a r \ ( v V a e ^ ^ 
+ — A f _ f (p)+ It-) A dust (p) 

+ AME(y) , 

where M and D represent monopole and dipolc residuals, 
respectively, A4 is the amplitude in antenna temperature 
units of the i th foreground component at the reference 
frequency {i/q, /^o, Ao), and j3, a and e describe the spec- 
tral response of synchrotron, free-free and thermal dust 
emission, and an AME contribution has also been in- 
cluded. 

In principle, there is no reason why these parameters 
should be constant over the sky, so we should allow them 
to vary pixel by pixel and let Commander solve for the 
most likely value. In practice, we are limited by the 
number of frequency bands, five in the case of WMAP, 
from which we also wish to determine the CMB contri- 



bution. Therefore, there remain only 4 maps that can 
be used to infer the foreground emission. Addressing the 
full problem as posed in Eq. [1] in an independent and 
self-consistent way is therefore not possible. Since our 
present analysis is motivated by a desire to assess the 
presence of a hard synchrotron component in the Galac- 
tic center, the microwave haze, we prefer to disentangle 
the foreground contributions in the low frequency range. 

A plausible and widely used alternative, which allows 
for more realistic foreground modeling, is the use of ex- 
ternal templates. At very low frequency (< 1 GHz), 
the observed sky signal is dominated by the synchrotron 
emission from our Galaxy, and is little contaminated by 
free-free emission, at least away from the Galactic plane. 
A full-sky map of the synchrotron emission is provided by 
the 408 MHz map of lHaslam et all (|1982fl . The optical 
Ha line is known to be a tracer of the free- free con t inuum 
emission at microwave wavelengths. iFinkbeinerl ()2003l ) 
produced a full-sky Ha map as a composite of various 
surveys in both the north e rn an d southern hemispheres. 
Finally, iFinkbeiner et al.l (|1999l ) predicted the thermal 
dust contribution at microwave frequencies from a se- 
ries of models based on the CO BE-DIRBE 100 and 240 
/xm maps tied to COB-E-FIRAS spectral data. We use 
predicted emission at 94 GHz from the preferred Model 
8 (FDS8) as our reference template for dust emission. 
These templates, smoothed to 60 arcminutes and down- 
graded to A s ido = 128, are shown in Figure [T] 

Since the thermal dust contribution to the WMAP 
bands is small compared to synchrotron and free-free 
emission, we can describe it by means of the FDS tem- 
plate with a fixed spectral index j3 = 1.7, allowing for 
an overall amplitude, b. The same can be done for the 
bremsstrahlung emission, assuming the Ha template as 
a sufficiently accurate description at 23 GHz and rescal- 
ing the amplitude according to a power law with index 
a = —2.15. The synchrotron emission is expected to vary 
across the sky, and the template we have is at 408 MHz, 
quite a large stretch from the first band of WMAP. A 
wise choice is to let Commander solve for an amplitude 
and spectral index at every pixel, choosing 23 GHz as 
pivot. It should be noted that variations in the gas tem- 
perature will be imprinted onto the synchrotron compo- 
nent, as will deviations in the spectrum of thermal dust 
(which is a much smaller contribution). 

A detailed study of the foreground emission in the fre- 
quency range 2 3-94 GHz emplo y ing e xternal templates 
can be found in iDickinson et al.1 ()2009l ) . Notice that the 
value of the spectral index of the dust, 1.7, is somewhat 
dependent on the specific method used to derive it. A 
direct fit to the predicted dust templates at WMAP fre- 
quencies using the FDS8 model yields a lower value of 
/? = 1.55, which is in agreement with what we find when 
applying template fitting procedure to WMAP maps (see 
Figure©. 

3. WMAP ANALYSIS 

Previous work on the Galactic Haze relied heavily on 
template regression techniques, raising a debate as to 
whether the assumption of constant spectral behavior 
across the sky is the cause of the excess emission near 
the Galactic Center. To address this issue, we follow an 
alternate approach and solve instead for a simple model 
describing Galactic emission with two power laws. One 
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Fig. 1. — Top to bottom, the three templates used to trace 
synchrotron (Haslam 408 MHz, top panel), free-free (Ha, middle 
panel) and dust, both thermal and spinning components (FDS, 
lower panel). 



is a low-frequency component with a falling spectrum to 
account for synchrotron and free-free emission, as well 
as AME. The second is a higher frequency component 
with a rising spectrum to represent thermal dust emis- 
sion. This appears to be well-motivated in studies of the 
WMAP Maxim um Entropy Meth od (MEM) foreground 
solutions, as in IPark et al.l (l200l . While we can solve 
for the amplitude and spectral index at every pixel at 
of the low-frequency component, where the variability of 
the signal is large, we fix the dust emissivity at high fre- 
quency to e = 1.7, consistent with previous applications 
of the Gibbs sampling technique to the WMAP data. 
We choose 22.8 GHz and 94 GHz as pivot frequency for 
the low-frequency and thermal dust component respec- 
tively. The drawback of this approach is that the three 
low-frequency emission mechanisms are combined into a 
single power law component that may not provide an 
fully adequate description of their physical complexity. 
However, we will attempt to separate these components 
post-sampling in Section 13.11 

We str ess that our approach dif fers from that fol- 
lowed by iCumberbatch et al.l (|2009f ). where Commander 
was only used to extract a CMB map while modeling the 



foreground components with templates. The main focus 
of that work was comparison between the WMAP In- 
ternal Linear Combination (ILC) CMB map, previously 
adopted in studies of the Haze, and the posterior mean 
CMB map derived from Gibbs sampling. The presence of 
the Haze was found to be stable with respect to the CMB 
map choice. In the present paper, we take full advantage 
of Commander to jointly sample from foregrounds and 
CMB. In addition, we choose a less aggressive smooth- 
ing of the WMAP maps and mask, both factors that are 
likely to enhance the signal-to-noise ratio of the Haze 
component. 

The Commander result for the CMB map is shown in 
Figure [2j As a diagnostic, we show the mean x 2 map 
in the lower panel of Figure [2j The number of pixels 
with x 2 > 20, corresponding to the 99.9% confidence 
level for 5 degrees of freedom, is 0.7%. Most of these 
pixels lie near point sources and we conclude that they 
correspond to point source bleeding outside the mask 
after convolution with a large beam. The overall lack of 
features in the \ 2 map is an indication of the goodness-of- 
fit of the model and that the residuals at each frequency 
are compatible with our noise description. 

We checked that the Commander power spectrum is con- 
sistent witli_t42e_l2estestimate provided by the WMAP 
team ([Larson et al.|[20TlT) at the la level up to i = 200. 
Beyond that, the Commander solution has significantly 
more scatter due to the regularizing noise added to the 
smoothed, binned maps (see Sec. l2.1[) . In addition, while 
the WMAP team also used Gibbs sampling techniques 
for low multiples, their high I power spectrum was ob- 
tained by applying a quadratic estimator to the cross- 
spectra of V and W bands. We emphasize that the 
agreement between the two power spectra below I = 200 
is what is important for our purposes since we are con- 
centrating on features in the map that are much larger 
than 1°. The possible presence of a residual monopolc 
and dipolc in the AT data has been taken into account 
and the mean values obtained are displayed in Table 
[TJ They are smal l and compatible with those found 
by iDickinson et~aTl (|2009f ) in the WMAP 5 -year data, 
though it is important to keep in mind that monopolc 
and dipole features become strongly co upled to fore- 
ground s in CMB analyses, as discussed bv lEriksen et al.l 
(|2008bD . 

TABLE 1 

Mean values (fJ.K) for the monopole and dipole residuals 

AT EVERY FREQUENCY. 



K- band Ka-band Q-band V-band W-band 

23 GHz 33 GHz 41 GHz 61 GHz 94 GHz 

"M 2.4 ±0.9 5.6 ± 1.0 3.5 ± 1.1 2.7 ± 1.0 3.9 EE 1.0 

D x -2.6 ±1.2 -3.4 ±1.2 -2.7 ± 1.2 -2.7 ±1.2 -3.0 ±1.2 

D y -3.6 ±0.9 -4.3 ±0.9 -4.4 ± 0.9 -3.5 ± 0.9 -4.0 ± 0.9 

D z 1.9 ±0.2 1.9 ±0.2 2.4 ±0.2 1.8 ±0.2 2.0 ± 0.2 



Foreground maps (left column) and associated errors 
(right column) for low frequency amplitude and spectral 
index measured at 23 GHz and the thermal dust contri- 
bution at 94 GHz are shown in Figure [3J 

It is clear that our derived low frequency component 
represents a combination of multiple emission mecha- 
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Fig. 2. — CMB Commander posterior average, rms and mean \ 2 
map. No particular features are present, meaning that the model 
is a very good fit to the WMAP 7-year data. The lower limit 
results 1.9, whereas the upper limit, 20, corresponds to 0.001 
probability for a \ 2 distribution with 5 degrees of freedom. Within 
the Gibbs sampling framework, we do not search for the maximum- 
likelihood (ML) point, but rather sample from the posterior: in 
practice, we add a scatter term around the the ML point, drawn 
from the noise model of the data. This results into a number that 
is distributed with 5 degrees of freedom. 



nisms, which we now attempt to disentangle. The low 
frequency component shows a strong correlation with 
thermal dust emission at 94 GHz as modeled by FDS, 
correlation which is interpreted as signature of a spin- 
ning dust contribution. The spectral index map, second 
row in FigurcO is particularly informative in this respect 
since departures from the prior assumed in Commander, 
— 3.0±0.3 are driven by the data. The mean value at high 
Galactic latitudes turns out to be a bit higher than the 
prior (~ —2.9) as consequence of the superposition along 
the line of sight of multiple components. It is remarkable 
how very bright free-free regions show up clearly, requir- 



ing a spectral index close to —2.15, which saturates the 
scale in Figure [4] The visual correlation of these diffuse 
gas clouds with those present in the Ha map is striking. 
Spectral index values lower than the prior are visible in 
regions where both synchrotron and dust correlated emis- 
sion are present, and distinguishing between the two is 
difficult. 

3.1. Template regression 

Following a common approach to the problem (see 
for e xample iDobler fc~F inkbcincr 2008b; iBottino et al.l 
120101 and references therein), we regress our amplitude 
maps against the three templates shown in Figure [TJ 

Assuming the amplitude map to be a linear combina- 
tion of these templates, the coefficients A\ are computed 
from a \ minimization of the form: 



x 2 = EE( s "(p)-E^w) ( 2 ) 

v p i 



A" 



dA v t 



d 2 X 2 
d[A»] 



A A- 



where Tj are the foreground templates and is the 
Commander amplitude solution. The explicit expression 
for the regression solution is given by: 



V 



N$(p) ' 



(3) 



where Ty and £?J describe the noise weighted correla- 
tion between foreground templates and frequency maps 
and templates, respectively. The error AA^ 1 is given 

See also (Fernan dez-Cerezo et al.l 12006: 



by V(2TM) 

IHildebrandt et al.|[2007l ). Since templates are not reli- 
able in the Gala ctic plane, the fit is performed outside 
the Kq85 mask (|Larson etalll20Tl . which removes the 
plane and detected point sources. In addition, the mask 
covers regions of high dust column density where extinc- 
tion makes our Ha template a poor tracer of free-free 
template emission. 

The coefficients of the fit are quoted in Table [2] and 
the resulting template and residuals in Figures [5] and HI 
Contrary to the usual technique, where the regression is 
performed on the WMAP maps, we regress out instead 
the Commander posterior average maps, which have noise 
contributed by the sampling procedure. Notice that the 
error we quote is obtained by integrating the posterior 
distribution over the other parameters, indices and com- 
ponent maps, and hence accounts for the uncertainty as- 
sociated with the foreground fit as well as the instru- 
mental noise. We also remind the reader that the input 



6 



Low-frequency Comp. 9 33 GHz: WMAP7 w- frequency Comp. Commander Posterior RMS: WMA! 
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Fig. 3. — Mean field map (left column) and rms (right column) for low-frequency component amplitude (upper panel) and spectral index 
(middle panel) and dust amplitude (lower panel) for the foreground model with two power laws. 
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Low-frequency Comp. Index @ 23 GHz: WMAP7 
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Fig. 4. — Full sky spectral index mean field map. The Galactic 
plane shows strong variation corresponding to regions where one 
component among synchrotron, free-free and spinning dust emis- 
sion dominates. 



maps analyzed by Commander are noisier than the ex- 
pected from simply smoothing them because we add reg- 
ularizing noise required by the sampling algorithm (see 
Section I2TT]) , These two factors increase the uncertainty 
on the template amplitudes. A detailed investigation is 
presented in Appendix |B| 

As expected, the high frequency foreground component 
correlates with the FDS model of thermal dust only. The 
low-frequency component is mainly a combination of syn- 
chrotron emission and correlated dust, the latter possi- 
bly due to spinning dust. We note that the contribution 
of free-free emission is weak because the Galactic plane, 
where it is strongest, has been masked out. To further 
quantify the goodness-of-fit, we show scatter plots of the 
mean field Commander solution as function of the derived 
linear template shown in Figure [7J 



TABLE 2 

Regression coefficients of the Commander foreground 

AMPLITUDE MAPS. WHILE THE THERMAL DUST CORRELATES WITH 

THE FDS MODEL ONLY, THE LOW-FREQUENCY COMPONENT IS 
INDEED A MIXTURE OF ALL THREE TEMPLATES, SYNCHROTRON AND 
DUST EMISSION BEING THE STRONGEST. MOREOVER, WHILE THE 
THERMAL DUST RESIDUALS ARE FEATURELESS AND CONSISTENT 
WITH NOISE, THE LOW-FREQUENCY COMPONENT RESIDUALS SHOW A 
CLEAR EXCESS OF POWER AROUND THE GALACTIC CENTER. 
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Fig. 5. — Low-frequency foreground amplitude maps (top) com- 
pared to the linear combination of templates (bottom). The dif- 
ference (i. e. residuals) is shown in the middle panel. 



WMAP frequencies outside the Galactic mask, 94 GHz 
being only slightly contaminated, the latter is poorly 
characterized and traced through the correlation between 



WMAP 7-yr : 


Dataset 


Haslam 


Ha 


FDS 


— 


Low Freq. Comp. 


(3.6 ±1.2) x 10" b 


6±6 


7±2 


0.94 I 


Thermal Dust 


(0.02 ± 1.0) x 10" 6 


-0.6 ±5 


1.1 ± 1.8 


0.58 C 



Although the correlation is clearly present and follows 
the y = x line, looking at the residual map is very in- 
structive. In the case of dust, the residuals arc compat- 
ible with noise and this explains why the points cluster 
close to y — line. The situation is more intriguing for 
the low frequency component where the fit is not per- 
fect and leaves an excess of power in the proximity of 
the Galactic center - the Galactic Haze - which has been 
advocated to be a distinct contribution. Fainter positive 
and negative regions exist as well, perhaps suggesting an 
overly-simplistic modeling of the dust component. 

4. FURTHER ANALYSIS 

The most troublesome foreground is the dust, both 
thermal and spinning. While the former remains week at 



0.58 our foreground model by first regressing WMAP chan- 
nels against Haslam, Ha and FDS to obtain a spectral 
energy distribution for each foreground, and then use the 
resulting phenomenological dust SED with Commanderto 
solve for the dust amplitude and a low-frequency compo- 
nent, described by a power law, together with the CMB. 
We recognize that this spinning dust model is simple and 
unlikely to capture the full complexity of the dust emis- 
sion. However, more complicated models, like spatially 
varying spinning dust index or multiple component spin- 
ning dust, although physically motivated (such as those 
proposed by ??), would require larger frequency coverage 
for practical use. We therefore adopt our simple model 
as an adequate description for the available WMAP data 
set. 



4.1. Low-frequency Component and Phenomenological 

Dust 
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Fig. 6. — Thermal dust amplitude maps (top) compared to the 
linear combination of templates (bottom) the difference (i. e. resid- 
uals) is shown in the middle panel. 



The first Commander run with this model returned a 
very steep spectral index for the low-frequency compo- 
nent, close to the limit of our prior for some regions in the 
Galactic plane known to have strong spinning dust emis- 
sion. The x 2 °f our solution in those regions was high, 
a sign that the foreground model fails. We then tuned 
the amplitude of our dust model, increasing it at low fre- 
quencies until we obtained a shallower spectral index for 
the low frequency component and smaller x 2 values for 
the same regions. In practice, this means that we mod- 
ify the SED template to increase the relative amount of 
spinning dust to thermal dust. The resulting spectral 
response is shown in Figure [8j together with those of the 
soft synchrotron (/? = —3), free-free (/? = —2.15) and 
thermal dust components (J3 = 1.7). 

Figure [5] shows the mean field amplitude maps of the 
thermal/spinning dust, and of the low-frequency com- 
ponent along with its spectral index. At high Galac- 
tic latitudes, where the signal-to-noise is very low, the 
spectral index is consistent with the Gaussian prior, 
j3 = — 3.0±0.3, —2.9 being the mean value, whereas closer 
to the plane the solution is driven by the data. In high 
latitude regions of strong synchrotron features, such as 
Loop I, the spectral index is noticeably softer than close 
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Dust @ 94 GHz: WMAP7 
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Fig. 7. — For the low- frequency foreground (top) and dust ampli- 
tude (bottom) we show the pixel-to-pixel comparison between the 
best-fitting template combination (y-axis) and the Commander val- 
ues (x-axis). The red line marks the y=x scaling. Blue and green 
dots denote residuals with a significance larger than 3<r, positive 
and negative, respectively. Orange points are pixels corresponding 
to the Haze: 3a positive excess within a disc of 36 degrees around 
the Galactic center. 
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Fig. 8. — SED for the phcnomcnological dust, compared to those 
of other components of interest: synchrotron, free-free and thermal 
dust emission. 



to the plane, where spinning dust and free-free emission 
become more important. Finally, it is interesting to note 
that the dust map recovered at 94 GHz is remarkably 
similar to the FDS prediction. 

As noted above, the low-frequency Commander solu- 
tion represents the combination of several different emis- 
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Low- frequency Comp. <8> 33 GHz: WMAP7 
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FlG. 9. — Low-frequency component mean field amplitude (upper 
panel) and spectral index (middle panel) map at 23 GHz and dust- 
correlated emission (lower panel). 



sion mechanisms known to coexist at 23-41 GHz: syn- 
chrotron, free-free, and possible spinning dust residu- 
als. This is confirmed by the low-frequency spectral in- 
dex map that presents shallower values at high Galactic 
latitude, but clearly emphasizes regions in the Galac- 
tic plane with a very steep spectral index, likely to be 
spinning dust clouds. A comparison between the low- 
frequency component spectral index and amplitude for 
the two models is shown in Figure [TU1 The evidence for 
dust-correlated emission is compelling, implying that our 
phcnomcnological dust SED is a good approximation of 
the sky. 

For the most part, all of these emission mecha- 
nisms decrease in intensity above 23 GHz (but see 
iDobler fc Finkbeinerl [2008bT iDobler eTail 12001) . and it 
is because of the noise in the data that a single power 
law results in a good \ 2 ■ ^ n an effort to disentangle these 
primary sources of emission, we again apply template re- 
gression to the Commander outputs using the Haslam 408 
MHz map, the Ha map, and the FDS map as tracers of 
synchrotron, free-free, and dust (thermal and spinning) 
emission, respectively. 

The regression coefficients arc quoted in Table [3] and 
indeed they confirm our success: outside the applied 
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Fig. 10. — Difference between two Commander posterior mean low- 
frequency component amplitudes (upper panel) and spectral in- 
dices (lower panel) assuming two different foreground models: two 
power-laws, and a power law with the phcnomcnological dust SED. 



mask, dust is positively correlated with the FDS map 
only, whereas the low-frequency component map can be 
described as a linear combination of Haslam and Ha. 
The goodness of the regression is expressed by the coef- 
ficient. 

TABLE 3 

Regression coefficients for Commander foreground 
AMPLITUDES: Top: WMAP 7-YR RUN; Bottom WMAP 7-YR AND 
Haslam. The Galactic emission is decomposed into physical 
components, each of which shows a clear correlation 
with one template only. 



Dataset 


Haslam 


Ha 


FDS 




WMAP 7- 


■yr 




Low-frcq. Comp. 


(3.4 ±3.3) x icr b 


9 ± 15 


0.6 ±6 


Dust 


(0.1 ± 1.3) x 10" 6 


-1±6 


2.3 ± 2.3 


Haslam+WJVLAP 7-yr 


Soft Sync. 


(5.9 ± 0.4) x 10" B 


-0.4 ± 1.3 


(-0.6 ±7) x lO" 1 


Low-freq. Comp. 


(-1.7 ±3.5) x 10- fi 


8 ± 15 


0.4 ±6 


Dust 


(-0.3 ± 1.2) x 10- 6 


-0.25 ± 6 


2.4 ±2.2 



As in the former run with thermal dust described by 
a power law, the FDS template shows itself a remark- 
able tracer of the dust map recovered by Commander, 
with the residual showing no evidence for large scale de- 
viations from the template. The same is true for the 
low frequency component with the significant exception 
of the "haze region" around the Galactic center. Fig- 
ure [11] shows the residuals after removing the template- 
correlated emission from the Commander mean field low- 
frequency foreground map. In other words, synchrotron 
emission elsewhere (in Loop I and in the Galactic plane) 
as well as free-free emission are traced well by our tem- 
plates while the haze is definitely not, indicating that this 



10 



Low- frequency Conip. <8> 33 GHz: WMAP7 




-200 400 ^.K Antenna 



iidual Foregrounds Hard Synchrotron II 23 GHz: WM^ 




mm Antenna 



Template Linear Combination: WMAPT 




Fig. 11. — WMAP 7-yr: Comparison between low-frequency fore- 
ground component amplitude solution and linear fit of the tem- 
plates: an excess of power around the Galactic center is present. 



large structure is not morphologically correlated with 
any of the templates. We find that the haze emission 
is present in the data and contributes significantly to the 
total emission towards the Galactic center (and particu- 
larly in the south). 

Our findings are in agreement with what has been 
described by severa l autho rs and in particular by 
iDobler &: Finkbeinerl (|2008al ). who claim the presence 
of an additional foreground component characterized by 
a spectral behaviour harder than a synchrotron compo- 
nent, but compatible with neither free-free (because the 
spectrum is too soft) nor spinning dust emission (because 
of lack of a thermal feature). Previous analyses were 
based on presubtraction of CMB cleaned maps from the 
data, which is problematic due to the fact that no CMB 
estimator is completely clean of foregrounds. This leads 
to a bias in the inferred spectrum of the Galactic emis- 
sions. To the extent that it is possible with five bands, 
we have attempted to reduce this systematic by simulta- 
neously solving for the cleaned CMB map and the fore- 
ground maps while also determining a spectral index on 
the sky, within a Bayesian framework, where the good- 
ness of the fit is controlled pixel-by-pixel through a \ 2 
evaluation. 



4.2. WMAP 7-year Data combined with Ancillary Data 

Our results in the previous section indicate that the 
haze represents either a new component not present in 
the Haslam map or a variation of the spectral index as- 
sociated with Haslam 408 MHz data. In order to as- 
sess this, we include the 408 MHz Haslam map together 
with WMAP data and run Commander on six bands cov- 
ering a factor of ^250 in frequency. This allows us to 
fit for an additional component and separate soft syn- 
chrotron from other emission mechanisms, the former 
being described by a power law with fixed spectral index 
p = —3, the latter by a spatially varying spectral index. 
We expect the soft contribution to be mainly driven by 
the Haslam map. Our foreground model in this second 
Commander run becomes: CMB, dust described by a spa- 
tially constant spectrum which accounts for both spin- 
ning and thermal emission, soft synchrotron with fixed 
spectral response and an additional low frequency com- 
ponent with spatially varying spectral index: 

T„(p) =M„+ D »AP) + ( — Y {p)Aowhe q (p) 

d—x,y,z 

+ {j^) Aoftsynch(p) + s{v)A dnst (jp), 

(4) 

where v Q = 22.8 GHz and A = 408 MHz, and s(v) is 
the fixed effective dust spectrum, normalized to /io = 
33 GHz. 

The mean posterior CMB we obtain with this 
Commander model is slightly different compared to the 5- 
band case showing less power close to the Galactic plane. 
This suggests that the improved flexibility of our model 
has resulted in better component separation. The over- 
all x 2 of these 6-band run does not change, since we add 
one map and we fit for one more foreground component. 
It is interesting to notice that we retrieve a value of the 
monopole in the Haslam map of ~ 3.2K. 

In Figure [T2] we report the mean field map of the 
foreground amplitudes: by visual inspection, the three 
amplitude maps look strongly correlated with the fore- 
ground templates. The coefficients of the regression are 
summarized in Table [3] 

The residuals obtained after regressing out the tem- 
plates from the Commander solution obtained for this 
more flexible foreground model are shown in Figure 1131 
The residuals arc quite striking and demonstrate that the 
soft synchrotron component (i.e., emission with a spec- 
tral index of t/~ 3 ) is almost perfectly correlated with 
Haslam and is not present, while the low frequency com- 
ponent clearly has a v~ 215 free-free component that is 
strongly correlated with Ha and the haze component, 
which is softer. This may suggest that the haze is not 
merely a spectral index variation of emission present in 
the Haslam map. In fact, the low frequency component 
is anti-correlated with Haslam (driven primarily by the 
small negative residuals of Loop I). It is plausible that 
the prior we set for the spectral index of the soft com- 
ponent, ft = —3, is not steep enough to characterize the 
Spur (J3 = —3.05), although at high Galactic latitudes 
a flatter behavior is measured (/3 ~ —2.9). Driven by 
this consideration, we ran Commander setting the spec- 
tral index of the synchrotron component to (3 = —2.9 
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Fig. 12. — WMAP 7-yr+Haslam: Commander posterior mean fore- 
ground amplitude maps. 



and P = —3.05: the former resulted into more nega- 
tive residuals in the synchrotron component, the latter 
reduced without completely removing them. 

The haze is a relatively weak residual compared to the 
soft synchrotron component , thus a spectral index ch ange 
is required to generate it ([Mertsch fc Sarkail I201C0 . If 
this were interpreted as a simple variation of the spectral 
index of Haslam rather than a separate component not 
visible in Haslam, then this would imply that the spec- 
tral index is harder from 23-60 GHz than from 0.408-23 
GHz. i.e., the spectrum would be concave up. However, 
cooling of electrons always generates spectra which are 
concave down, because energy losses of cosmic-ray elec- 
trons are proportional to the electron energy squared. 
Thus, power-law injection of electrons results in concave 
down spectra as they cool. In other words, a different 
population of electrons with a harder spectral index must 
be present: a simple variation of the spectral index of the 
same electron population is not what we measured, since 
the amplitude of the soft component that we retrieved is 
strongly correlated with the Haslam map at 408 MHz. It 
is true that the assumption of a power-law for the syn- 
chrotron component can be too simplistic, and a (nega- 
tive) curvature might be present, as a result of a break 
in the spectral behaviour. However, given the number of 



Fig. 13. — WMAP 7-yr+Haslam: Comparison between low- 
frequency component amplitude solution and linear fit of the tem- 
plates. The Haslam-uncorrclatcd low frequency component shows 
an excess of power which is not present in any of the templates. 



frequencies available, sparse between 0.408 and 23 GHz 
and relatively narrow range of the WMAP data, allowing 
an extra parameter to be determined is not possible. 

Again, the dust component is only strongly correlated 
with the FDS map, even more so than in the previous 
model, again likely indicating a better foreground sepa- 
ration. 

Lastly, we point out that the spectral index map of 
the haze clearly shows that its spectral behaviour is not 
compatible with free-free, but it is driven by the prior 
(/3 = —2.5 ± 0.3). Unfortunately, this may suggest that 
the regularizing noise we added is too high to enable us 
to be definitive about the spectral index of the haze. 

Up until now we have explored template fits in which 
the spectra are taken to be uniform across the sky. 
This is obviously not the case, and so to demonstrate 
the effects of this assumption we break the sky into 
regions and perform individual fits on these. Exam- 
ples o f t his technique can be f o und in iHildebrandt et al.l 
(120071): IDobler fc FinkbeTnen (I2008aj) : IDoblerl (1201 HI : 
iGhosh et al.l ( 2012I ). In practice, it means computing 
the correlation coefficients given in Eq. Hon a subset of 
pixels for the chosen regions. Here we will not perform 
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such a detailed analysis, but restrict ourselves only to 
the Haze region around the Galactic center. The idea 
behind this approach is that if the Haze corresponds to a 
region of global variation of the spectral index of the syn- 
chrotron component, the correlation between our resid- 
ual map and the Haslam map computed on the Haze 
pixels would be high. We can select the brightest pixels 
based on the rms map we derived from Commander: in 
Figure [14] we show the significance of the residuals and 
the pixels above 3cr which we chose within 36 degrees of 
the Galactic center. The correlation coefficients we ob- 

al Foregrounds Low-freq, Comp. @ 23 GHz; Haslam+ 
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Galactic Haze Region at 23 GHz 
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Fig. 14. — Significance of the residuals for the soft synchrotron 
component (top panel) and pixel selected for the study of the cor- 
relation on a small region around the Galactic center. 



tained for the Haze are 3.3 ± 3.5 x 1CT 6 , 0.1 ± 21 and 
3 ± 7 respectively for Haslam, Ha and FDS templates, 
showing no evidence for strong correlation. It is interest- 
ing to notice that while the amplitudes for free-free and 
dust emission remain stable with the threshold, this is 
not the case of the Haslam map, whose correlation de- 
creases with the lower bound of the chosen significance 
cut. We tried 2.0, 2.5, 3.0, 3.5 a finding more and more 
correlation: this is not surprising since the number of se- 
lected pixels decreases and the statistics degrades. This 
is a due to the regularizing noise. 

This approach will be more effective when larger fre- 
quency coverage becomes be available. 

To improve the constraints on the spectral index of 
the hard low-frequency component, we re-ran comman- 
der without adding noise to the input maps and adopt- 
ing a different prior, f3 = —2.15 ± 0.3, with the idea that 
the higher signal-to-noise ratio would allow the data to 
drive the solution. (In order to speed converge, we also 
dropped the sample of C;s). Except for regions with 
strong free-free emission, we now observe a mean spec- 
tral index of — 2.3±0.27 for |6oai| > 30°, whereas we find 
—2.4 ± 0.22 for the haze region. Though not conclusive, 



this docs support the hypothesis of a har der spectral in - 
dex for the haze component proposed by iDoblerl (|2011l) . 
The regression coefficients for this case are 0.4±2 x 10~ 7 , 
9 ± 7 and 0.08 ± 2 for the Haslam, Ha and FDS tem- 
plate respectively, which show a clear improvement of the 
component identification. Figure [15] shows the posterior 
mean amplitude and spectral index of the low-frequency 
component with a 35° latitude around the Galactic cen- 
ter. 



Haze Region: Amplitude 
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Fig. 15. — Haze amplitude and spectral index obtained after 
regressing Commander posterior mean amplitudes found without 
adding rcgularization noise. 



5. CONCLUSIONS 

We applied the Gibbs sampling technique implemented 
in the Commander software to WMAP 7-year data aim- 
ing at characterizing, simultaneously, CMB (map and 
angular power spectrum) and foreground emission. Our 
analysis improves previous work in a number of aspects. 
We pushed the analysis to higher angular resolution, 1°, 
which enabled us to directly compare the foreground 
solutions to available templates, resulting in excellent 
agreement. This represents a success on its own but also 
confirms the power of the Bayesian approach, which is 
likely to perform better when the larger Planck data set 
becomes available. 

Higher angular resolution requires better modeling of 
the noise in the maps after the smoothing procedure: we 
added scan-modulated regularizing noise estimated via 
Monte Carlo simulations. Despite the limited number of 
frequencies available, we did not directly use foreground 
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templates in the Gibbs sampling run, but rather let the 
parameters of the model vary across the sky. This al- 
lowed us to distinguish regions in the sky where a specific 
foreground mechanism dominates on the basis of the its 
spectral behaviour. The presence of a strongly dust cor- 
related emission at low frequencies, explained by invok- 
ing spinning dust grains, emerges not only through the 
amplitude map but also in the spatial variation of the 
spectral index of the low frequency component, which 
mainly results from a mixture of synchrotron emission 
and dust-correlated emission. The number of available 
frequencies still limits us, since it forces us to use a con- 
stant spectral index for the dust. 

Regressing Commander solutions against foreground 
templates is a complementary way to disentangle 
the various emission mechanisms. It is by apply- 
ing this procedure that we confirm the presence of 
excess of signal localized around the Galactic cen- 
ter, the microwave haze. A hint of s uch a emis- 
sion has b een found in the Fermi data (|Dobler et al.l 
120101 I2011D . The study of this component has stim- 
ulated a rich production, both pointing out possible 
systematics in the applied regression procedures (,see 



for in stance ILinden fc Profumo 120101 : IMertsch fe Sarkarl 
l2OT0l and) and proposing possible physical mecha- 
nism s (TBottino et al.ll2010l:lDobler et alil201lt iGuo et al.l 

TIT 
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The WMAP team remains ske ptical about the presence 
of the haze (|Gold et al.l 120111) . their main counter ar- 
gument being the lack of evide nce for such a signal 
in the polarization maps (though iDoblerl 120111 showed 
that such a signal would likely not be detectable in 
the WMAP data given the noise). Regarding the spa- 
tial correspondence between the WMAP haze and the 
Fermi haze/bubbles, we note that the increased noise 
in our output maps makes a direct comparison difficult. 
Comparison with e a rlier WMAP da ta was presented by 
iDobler et all (|2010l) ; ISu et al l (l2010f) and more recently 
with WMAP 7-year data bv IDoblerl (poTTh . We find, as 
in those studies, that the morphological correlation be- 
tween the two is reasonable at |6cai < 30°, where the 
signal in microwaves is most unambiguous. 

We emphasize that our analysis overcomes the problem 
of CMB subtraction and the circularity argument aris- 
ing from the removal of an internal linear combination 
of the channels used in previous analyses. Moreover, we 
addressed the issue of possible contamination of spinning 
dust as an explanation of the excess signal by including a 
correlation between the low frequency emission and ther- 
mal dust. We computed the spectral energy distribution 
of the dust as described by the FDS model, performing a 
template fitting of the WMAP channels as a preprocess- 
ing step. We input the resulting SED to Commander, solv- 
ing for an amplitude map, together with a low-frequency 



component. The resulting dust map is highly correlated 
with the FDS model, and the low-frequency component 
results from a sum of synchrotron and free-free emission 
only. When regressing this map against foreground tem- 
plates, the Galactic Haze is still present. One criticism 
to our approach could be the simplicity of the spinning 
dust model, but since one of our goals was to be indepen- 
dent of external data sets within a completely Bayesian 
framework, adding more complexity to the dust model is 
not possible given the available degrees of freedom. 

We also argued that the haze is not an artifact of the 
spectral variation of the synchrotron component across 
the sky. To this end, we included the Haslam map at 408 
MHz in the Commander run and expanding the foreground 
model: dust, combing both thermal and spinning contri- 
bution, synchrotron emission with fixed spectral response 
j3 = —3 and an additional component with a free spec- 
tral index, which would describe free-free, and any other 
contribution. We found that we could completely sepa- 
rate synchrotron and dust emission, and that we were left 
with an amplitude map that can be easily characterized 
by its spectral index. Free-free emission clearly shows 
up in the Galactic plane, together with the Haze which 
seems to have a harder spectrum than the synchrotron 
component. Unfortunately, outside the Galactic plane, 
the index map is quite noisy and the haze does not have 
a distinctive signature compared to the rest of the sky, al- 
though it is definitely not compatible with free- free emis- 
sion. 

Our analysis has improved previous knowledge of the 
foreground and shed new light on the nature of the Galac- 
tic Haze. We addressed at least three criticisms often 
attached to haze studies: coherent CMB removal, con- 
tamination from spinning dust, and the spatial variation 
of the synchrotron component spectral index. The ex- 
cess of signal seems to be stable with respect to them. 
However, we have made assumptions on the dust spectral 
behaviour, forced by the number of frequencies available, 
and we did not rule out the possibility of the presence of 
curvature in the synchrotron component spectral index 
or a model with a broken power law. Planck data will ad- 
dress these remaining issues and enable a more complete 
foreground model. 
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APPENDIX 

THE CMB MAP AND POWER SPECTRUM POSTERIOR 

Here we review the basic concept behind Gibbs sampling. Let us first focus on the case of one frequency map and 
no foregrounds. The data model for this case is 

d = s + n, (Al) 

where d is the data, s the CMB sky signal, and n instrumental noise. We assume both the CMB signal and noise to 
be Gaussian random fields with vanishing mean and covariance matrices S and N, respectively. The CMB sky can be 
written in spherical harmonics as s = m a tmXimi with the CMB covariance matrix then fully characterized by the 
angular power spectrum Ct according to Ce. m ,i'm' = {a*i m a-t'm') = Ci5u'S mm > . The noise matrix N is left unspecified 
for now, but we note that for white noise it is diagonal in pixel space, N%j = erf Si j, for pixels i and j and noise variance 

Our goal is to sample from the posterior density for both the sky signal s and the power spectrum Cg, given by 

P(s, C t \d) oc P(d|s, cyP(s, C t ) (A2) 
cc P{&\B,C t )P(*\Ct)P{Ct), (A3) 

In what follows we assume the prior P(Ct) is uniform. Since we have assumed Gaussianity, the joint posterior 
distribution may be written as 

21+1 °£ 

P(s, C t \d) c< e-iCd-^N-^d-.) Y[ e 2 2t+ T l P(C t ), (A4) 

t cv 

where we have defined the quantity at = -^I+i Yf m =-e \ a £m\ 2 as the angular power spectrum of the full-sky CMB 
signal. 

For the case here with the CMB signal assumed to be a Gaussian field, one can integrate over the CMB sky signal 
and analytically solve for the marginalized posterior P{C/\d). However, evaluating the posterior numerically for any 
specific angular power spectrum is computationally prohibitive as it involves the computation of the inverse and 
determinant of very large matrices. We therefore sample from the posterior using a Gibbs sampling algorithm. 
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Gibbs Sampling 

One procedure to sample from the joint density P(s, C^|d), as proposed by IJewell et al.l (|2004D and lWandelt et all 
(|2004l ). is to alternately sample from the respective conditional densities 

s l+1 <-P(s|C5,d) (A5) 

C\ +1 ^P(C e \s l+1 ,d). (A6) 

Here <— indicates sampling from the distribution on the right-hand side. After some "burn-in" period, the joint samples 
(s l ,CJ) will be distributed from the joint posterior. Thus, the problem is reduced to that of sampling from the two 
conditional densities P(s\d,d) and P(Ce\s,d). 

The conditional density P(Ce\s, d) in this case is independent of the data, P(Cg\s, d) = P(Ce\s), simply because the 
underlying CMB sky signal provides all the information needed to estimate the ensemble angular power spectrum Ci . 
Under the assumption of Gaussianity and isotropy, this conditional is given by the inverse Gamma distribution, Sf. 

1 to-l 2t + 1 Zi 

P(C e \s) <x =- = M+1 . (A7) 

In order to sample from this conditional density we first draw 2£ — 1 normal random variates pi, compute the sum 
p\ = Y?k=\ \Pe\ 2 > an d finally set 

Ci = (A8) 
PI 

giving a sample distributed according to the inverse Gamma distribution. 

The conditional density for the sky map given the angular power spectrum and data follows directly from the form 
of the joint Bayes posterior in Equation IA41 and given by 

P(s\Ci, d) oc P(d|s, C e )P(s\C £ ) (A9) 
K g-^d-^N-^d-.) e-^s-i* (A10) 

OC e -|(s-§) t (S- 1 +N- 1 )(s-s) ) (An) 

where we have defined the so-called mean-field map (or Wiener filtered data) s = (S _1 + N _1 ) _1 N _1 d. Thus, 
P(s\Ce, d) is a Gaussian distribution with mean equal to s and a covariance matrix equal to (S _1 + N -1 ) -1 . In order 
to sample from this conditional, we first generate two independent white noise maps luq and u>\, and solve 

[S^ 1 + INT 1 ] s = N _1 d + S^u; + N _ ^ Wl) (A12) 

The resulting map s is exactly a sample from the conditional P(s|Ce, d). The addition of the white noise maps simply 
reflects our uncertainty in the true but unobserved CMB sky - there are many CMB maps that are consistent with 
the data and power spectrum estimate, and we are simply making a random choice from this set of maps. 

The discussion so far was limited to a single band and no m odeling of instrumental response. The generalization 
is straightforward and can be found in lEriksen et al.1 (|2007bl ) , together with a discussion on the actu al num erical 
implementation which reduces the round-off errors. Here we quote the main result, which generalizes Eq. IA12l to: 

S-i + ^AtN^A, 

(A13) 

where A„ describes the beam response of the detector v and N its noise properties. Note that we now draw one white 
noise map for each frequency band, w v . The sampling procedure for P(C(\s) is unchanged. 

The foreground sampler 

The Gibbs sampling algorithm can be easily extended to include foregrounds described by a parametric model, F(Q). 
Including the foreground model parameters into the joint posterior, the sampling procedure generalizes to: 

s' i+1 «- P(s|Q,6>\d) 

C} +1 ^ P(C e \s l+1 ,d) 

^ p(0|C| +1 ,s i+1 ,d). 

We note that for the foreground models of interest, we do not have an algorithm to produce an exact sample from the 
conditional density P{9\ Cl +1 , s' +1 , d) but we i nstead generalize Gibbs sampling to an MCMC algorithm, the details 
of which can be found in lEriksen etaLl (j2008b[ )). 
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The parametric family of data models including foregrounds implemented in Commander are of the form 

M 




(A14) 



n=l k=l 



We may identify three main classes of foregrounds: 

i) t m are M templates multiplied by an amplitude at every frequency, a v ,m'i 

ii) f rl arc N templates whose spectral behaviour is known and described by the function f n {v)] we allow for an 
overall rescaling, b n , for each of them; 

iii) K foregrounds are described by a map of coefficients, Cfc, multiplied by the spectral response which is function 
of the frequency and the parameters of the foreground model, Ok- 

We remind the reader that the bold face notation means an array of size iVpi x . An example of the first class of 
foregrounds is given by monopole and dipole residuals; a special case of the second class is a free-free template 
(Figure [TJ top panel) whose spectral behaviour is known and follows the relation (v / i/q)~ 2A5 ; for the third type we 
may quote synchrotron emission described by an amplitude we solve for at the reference frequency, e. g. /io = 23 GHz, 
and a spectral response given by (v/hq)". 

It is instructive to compute the degrees of freedom for such a foreground model applied to WMAP data. We solve 
three maps, s, j3 and A sync h, 5 monopoles and dipoles, and two overall amplitudes if we describe free-free emission and 
thermal dust by means of the second class of foreground models. In total we look for 3 -/V p i x + 22 parameters. This is 
already pretty close to the maximum number of parameters allowed by the WMAP, < 5 iVpi x . We could ask for one 
more map, either dust or free-free, assuming a single power law and fixing the spectral index, or allowing a curvature 
term in the synchrotron model. Since a simple power law for thermal d ust and free-free emissi on has been shown to 
be consistent with the data in the frequency range spanned by WMAP ((Dickinson e t al. 200^), this may be used to 
disentangle the spinning dust contribution from thermal dust. 

We notice that to solve for the spectral index of a given component, the instrumental beam response of each channel 
must be taken into account. Up to now, Commander is able to work with maps at the same angular resolution only. 
Smoothing all fre quency maps and anci llary data to a common angular scale is then a necessary pre-processing step. 

As discussed in lEriksen et al.l (l2008bl ). it turns out to be more efficient to sample all the map amplitudes at once, 
followed by the spectral response parameter, and finally the angular power spectrum, following the iterative scheme: 



Sampling from the conditional density P(s, a m ,b n ,Ck\Cl,8 l ,d) is a generalization of Eq. IA131 whereas the sampling 
of the angular power spectrum remains unchanged, since Cg are functions of the sky signal only. Sampling of the non- 
linear degrees of freedom is through a standard inversion sampler: first compute the conditional probability density 
P(x\0), where x is the currently sampled parameter and denotes the set of all other parameters in the model, assuming 
the likelihood to be independent pixel by pixel: — 21n£(x) = X 2 = Y^ v {dv — s v {X, 9)) 2 /a 2 . Then, the corresponding 
cumulative distribution is computed, F{x\9) = J_ OQ P(y\6)dy, and the value of the x variable is chosen by drawing a 
uniformly distributed random number, u, and reading F(x\0) = u. 

This concludes the review on the implementation of Gibbs sampling as implemented in the computer code Commander. 



We have observed that the regression coefficients we found when fitting foreground templates to Commander posterior 
mean amplitudes are consistent with those discussed in other works, but they have much larger errors. We argued in 
the text that this is the result of the sampling procedure and has two causes: i) our input maps are noisier because 
of the additional noise term and ii) our uncertainties on the foreground amplitudes take into account the error on the 
other parameters of the model: CMB, other foreground amplitudes and spectral indices. In this respect, our errors 
are more conservative. 

To clearly show this, we perform the same regression described in Equations |3] and |3] on two different maps: 1) 
smoothed WMAP K-band and 2) Commander input K-band, which have to be compared to Commander output amplitude 
performance. The three maps are shown in Figure 1161 together with the fit residuals and the corresponding \ 2 ■ 

Table 0] compares the regression coefficients for the signal at 23 GHz. Consistently, the amplitudes are the same but 
the error bars increase dramatically due to the noise added to the maps. In particular we move from 10-30cr detection 
to l-3cr, which is driven by the scaling factor applied to the K-band (16, see Section |2T|) . 




l ,a m ,b n ,c k y +1 <- P(s,a m ,b n ,c k \C\,6\A) 
4- P(C e \s l+ \d). 



(A15) 



NOISE IMPACT ON THE REGRESSION PROCEDURE 
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This comparison suggests that a lower level of noise added to the input maps is useful to better characterize diffuse 
foreground emission. We will further investigate this issue in a forthcoming work. 

TABLE 4 

Regression coefficients of the Commander foreground amplitude maps compared to those obtained when regressing 

SMOOTHED WMAP CHANNEL AND Commander INPUT WMAP MAPS 



WMAP 7-yr 


Dataset 


Haslam 


Ha 


FDS 


r 


Smoothed K-band 
Commander Input K-band 
Commander output @ 23GHz 


(3.7 ±0.17) x 10~ b 
(3.7 ± 1.1) x 10- 6 
(3.6 ±1.2) x 10~ 6 


6.4 ± .6 
6±5 
6±6 


7.0 ± .25 
7±2 
7±2 


0.91 
0.85 
0.94 



Map @ 23 GHz Residual Foregrounds 




-200 400 fj,K CMB -200 200 pft Antenna 0.0 16.0 



FlG. 16. — Example of three regressions performed on the WMAP K-band smoothed to 60 arcminute (top), Commander input K-band 
(middle row) and output foreground amplitude (bottom). 



